Non-Integrable Galactic Dynamics 
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1 INTRODUCTION 

Galaxies have traditionally been viewed as integrable or nearly integrable systems, in 
which the majority of stellar orbits are regular, respecting as many integrals of motion as 
there are degrees of freedom. Three arguments have commonly been cited in support of 
this view. First, many reasonable potentials contain only modest numbers of stochastic 
orbits []. This is always true for the potentials of rotationally symmetric models, and there 
is even a class of non-axisymmetric potentials for which the motion is globally integrable, 
including the famous "perfect ellipsoid" (Kuzmin 1973; de Zeeuw & Lynden-Bell 1985). 
Second, stochastic orbits often behave in ways that are very similar to regular orbits 
over astronomically interesting time scales. Therefore (it is argued) one need not make a 
sharp distinction between regular and stochastic orbits when constructing an equilibrium 
model. Third, following the successful construction by Schwarzschild (1979, 1982) of 
^ ■ self-consistent triaxial equilibria, it has generally been assumed that the regular orbits - 
which are confined to narrow regions of phase space and thus have definite shapes - are 
the fundamental building blocks of real galaxies. 

Schwarzschild's discovery that many orbits in non-axisymmetric potentials are effec- 
tively regular came as a surprise, since triaxial potentials admit only one classical integral 
of the motion, the energy. In fact a modest fraction of the orbits in Schwarzschild's models 
were subsequently shown to be stochastic (Merritt 1980; Goodman & Schwarzschild 1981), 
though only weakly. But it was clear early on that certain modifications of Schwarzschild's 
potential could lead to a much larger fraction of chaotic orbits. For instance, Gerhard 
& Binney (1985) showed that the addition of a central density cusp or "black hole" (i.e. 
point mass) to an otherwise integrable triaxial model would render most of the center- 
filling, box orbits unstable, due to deflections that occur when a trajectory comes close to 
the center. This insight was followed by the discovery (Crane et al. 1993; Ferrarese et al. 
1994) that stellar spheroids generically contain power-law cusps in the luminosity density 

1 The terms "stochastic" and "chaotic" will be used interchangeably here. 



rather than constant-density cores. Evidence for central supermassive black holes also 
gradually accumulated (Kormendy & Richstone 1995). It is now believed - not only that 
black holes are universal components of galactic nuclei - but that their masses are pre- 
dictable with high precision given the global properties of their host spheroids (Ferrarese 
& Merritt 2000; Merritt & Ferrarese 2001). Thus it is no longer possible to discuss galaxy 
dynamics in terms of idealized models like Schwarzschild's with finite central densities. 

Non-integrability has two important consequences. First, some orbits in non-integrable 
potentials respect fewer isolating integrals than there are degrees of freedom. Such or- 
bits are chaotic and behave in ways that are very different from regular orbits: they are 
exponentially unstable to small perturbations, and occupy a phase-space region of larger 
dimensionality than the invariant tori of regular orbits. The time-averaged shape of a 
chaotic orbit is similar to that of an equipotential surface and hence such orbits are much 
less "useful" than regular orbits for reinforcing the shape of the galaxy's figure. Second, 
while regular orbits generally still exist in potentials that are not globally integrable, they 
are strongly influenced by resonances between the frequencies of motion in different di- 
rections. These resonances are present even in globally integrable potentials but have no 
effect on the structure of phase space; in non-integrable potentials, however, the reso- 
nances divide up phase space into alternating regions of regular and chaotic motion, with 
the lowest-order resonances "capturing" the largest parts of phase space. Most regular 
orbits in non-integrable potentials can be associated with a definite resonance and have 
a shape that reflects the order of the resonance. 

This article reviews the following topics: (1) Torus construction, a set of techniques 
for characterizing regular motion in non-integrable potentials and for detecting departures 
from integrability; (2) resonances and their effect on the structure of orbits; (3) the orbital 
content of triaxial potentials with central point masses; (4) mixing, the process by which 
the phase-space density of stellar systems approaches a steady state; and (5) the relation 
between chaos in the gravitational iV-body problem and chaos in smooth potentials. 



2 TORUS CONSTRUCTION 

In systems with a single degree of freedom, constancy of the energy allows the momentum 
variable p to be written in terms of the coordinate variable q as H (p, q) = E, and the 
dependence of both variables on time follows immediately from Hamilton's equations. In 
general systems with N > 2 degrees of freedom (DOF), such a solution is generally not 
possible unless the Hamilton- Jacobi equation is separable, in which case the separation 
constants are isolating integrals of the motion. An isolating integral is a conserved quan- 
tity that in some transformed coordinate system makes dH/dpi = f{qi), thus allowing 
the motion in q { to be reduced to quadratures. Each isolating integral restricts the di- 
mensionality of the phase space region accessible to an orbit by one; if there are N such 
integrals, the orbit moves in a phase space of dimension 2N — N = N, and the motion 
is regular. The iV-dimensional phase space region to which a regular orbit is confined 
is topologically a torus (Figure 1). Orbits in time- independent potentials may be either 
regular or chaotic; chaotic orbits respect a smaller number of integrals than N - typically 
only the energy integral E. Although chaotic orbits are not confined to tori, numerical 
integrations suggest that many chaotic trajectories are effectively regular, remaining con- 



fined for long periods of time to regions of phase space much more restricted than the full 
energy hypersurface. 




Figure 1. Invariant torus defining the motion of a regular orbit in a two-dimensional 
potential. The torus is determined by the values of the actions J\ and J 2 ; the position of 
the trajectory on the torus is defined by the angles d\ and 9 2 , which increase linearly with 
time, 9i = uj{t + 6®. 

The most compact representation of a regular orbit is in terms of the coordinates on 
the torus (Figure 1) - the action-angle variables (J,#). The process of determining the 
map (x, v) — » (J, 6) is referred to as torus construction. There are a number of contexts 
in which it is useful to know the (J, 9). One example is the response of orbits to slow 
changes in the potential, which leave the actions (J) unchanged. Another is the behavior 
of weakly chaotic orbits, which may be approximated as regular orbits that slowly diffuse 
from one torus to another. A third example is galaxy modeling, where regular orbits are 
most efficiently represented and stored via the coordinates that define their tori. 

Two general approaches to torus construction have been developed. Trajectory- 
following algorithms are based on the quasi-periodicity of regular motion: Fourier de- 
composition of the trajectory yields the fundamental frequencies on the torus as well as 
the spectral amplitudes, which allow immediate construction of the map — > x in the 
form of a Fourier series. Iterative approaches begin from some initial guess for x(#), which 
is then refined via Hamilton's equations with the requirement that the 0, increase linearly 
with time. The two approaches are often complementary, as discussed below. 



2.1 Regular Motion 

In certain special potentials, every orbit is regular; examples are the Kepler and Stackel 
potentials. Motion in such globally-integrable potentials can be expressed most simply 
by finding a canonical transformation to coordinates (p, q) for which the Hamiltonian is 
independent of q, H = H(j>); among all such coordinates, one particularly simple choice 
is the action-angle variables ( Jj, in terms of which the equations of motion are 

Ji = constant, 

OH 

6i = ujit + 01 ^ = — , i = l,...,N (1) 



(Landau & Lifshitz 1976; Goldstein 1980). The trajectory x(J,#) is periodic in each of 
the angle variables 0$, which may be restricted to the range < Q± < 2n. The Jj define 
the cross-sectional areas of the torus while the Qi define the position on the torus (Figure 
1). These tori are sometimes called "invariant" since a phase point that lies on a torus at 
any time will remain on it forever. 

Most potentials are not globally integrable, but regular orbits may still exist; indeed 
these are the orbits for which torus construction machinery is designed. One expects that 
for a regular orbit in a non- integrable potential, a canonical transformation (x, v) — > (J, 6) 
can be found such that 

j. = 0, 6i = u u i = l,...,N. (2) 

However there is no guarantee that the full Hamiltonian will be expressible as a continuous 
function of the Jj as in globally integrable potentials. In general, the map (x, v) — > (J, 9) 
w ill be different for each orbit and will not exist for those trajectories that do not respect 
N isolating integrals. 

The uniform translation of a regular orbit on its torus implies that the motion in any 
canonical coordinates (x, v) is quasi-periodic: 

x(i) = ^2 X fe( J) ex P \i> {IkVi + m k u} 2 + n k u 3 ) t] , 

k 

v(t) = ^2 v fe(J) exp [i (l k ui + m k u) 2 + n k u 3 ) t] , (3) 

k 

with (l k ,m k ,n k ) integers. The Fourier transform of x(t) or v(t) will therefore consist of a 
set of spikes at discrete frequencies u k = l k 0J\ + Tn k uo 2 + n k uo 3 that are linear combinations 
of the N fundamental frequencies u>i, with spectral amplitudes X fc (J) and V fc (J). 

2.2 Trajectory-Following Approaches 

The most straightforward, and probably the most robust, approach to torus construction is 
via Fourier analysis of the numerically-integrated trajectories (Percival 1974; Boozer 1982; 
Binney & Spergel 1982, 1984; Kuo-Petravic et al. 1983; Eaker et al. 1984; Martens & Ezra 
1985). The Fourier decomposition of a quasiperiodic orbit (equation yields a discrete 
frequency spectrum. The precise form of this spectrum depends on the coordinates in 
which the orbit is integrated, but certain of its properties are invariant, including the iV 
fundamental frequencies uii from which every line is made up, u k = l k u>i + m k u)2 + n k u 3 . 
Typically the strongest line in a spectrum lies at one of the fundamental frequencies; once 
the uJi have been identified, the integer vectors (l k ,m k ,n k ) corresponding to every line 
uj k are uniquely defined, to within computational uncertainties. Approximations to the 
actions may then be computed using Percival's (1974) formulae; e.g. the action associated 
with 9\ in a 3 DOF system is 

J 1 = ^l k (l k uj 1 + m k uj2 + n k uj 3 )\^ k \ 2 . (4) 

k 

Finally, the maps (6 — > x) are obtained by making the substitution uj{t — > 9i in the 
spectrum, e.g. 



x(t) = \ X k (J) exp [i (l k ui + m k u 2 + n k uj 3 ) t] 



= ^2 x k{J) exp [i {l k 9i + m k 9 2 + n k 9 3 )] 
k 

= x(e 1 ,e 2 ,e 3 ). 



(5) 



Trajectory-following algorithms are easily automated; for instance, integer programming 
may be used to recover the vectors (lk,rrik,nk) (Valluri & Merritt 1998). 

Binney & Spergel (1982) pioneered the use of trajectory- following algorithms for galac- 
tic potentials. They integrated orbits for a time T and computed discrete Fourier trans- 
forms, yielding spectra in which each frequency spike was represented by a peak with finite 
width ~ 7r/T centered on u>k- They then fitted these peaks to the expected functional 
form Xk sin[(cj — uJk)T]/(uj — ujk) using a least-squares algorithm. They were able to re- 
cover the fundamental frequencies in a 2 DOF potential with an accuracy of ~ 0.1% after 
~ 25 orbital periods. Binney & Spergel (1984) used equation (4) to construct the "action 
map" for orbits in a principal plane of the triaxial logarithmic potential. Carpintero & 
Aguilar (1998) have applied similar algorithms to motion in 2- and 3 DOF potentials. 

The accuracy of Fourier transform methods can be greatly improved by multiplying 
the time series with a windowing function before transforming. The result is a reduction 
in the amplitude of the side lobes of each frequency peak at the expense of a broadening of 
the peaks; the amplitude measurements are then effectively decoupled from any errors in 
the determination of the frequencies. Laskar (1988, 1990) developed this idea into a set of 
tools, the "numerical analysis of fundamental frequencies" (NAFF), which he applied to 
the analysis of weakly chaotic motion in the solar system. Laskar's algorithm recovers the 
fundamental frequencies with an error that falls off as T -4 (Laskar 1996), compared with 
~ T~ l in algorithms like Binney & Spergel's (1982). Even for modest integration times 
of ~ 10 2 orbital periods, the NAFF algorithm is able to recover fundamental frequencies 
with accuracies of ~ 10~ 8 or better in many potentials. The result is a very precise 
representation of the torus (Figure 2). 

Since Fourier techniques focus on the frequency domain, they are particularly well 
suited to identifying regions of phase space associated with resonances. Resonant tori 
are places where perturbation expansions of integrable systems break down, due to the 
"problem of small denominators". In perturbed (non-integrable) potentials, one expects 
stable resonant tori to generate regions of regular motion and unstable resonant tori to give 
rise to chaotic regions. Algorithms like NAFF allow one to construct a "frequency map" 
of the phase space: a plot of the ratios of the fundamental frequencies (uji/uj 3 , u 2 /uo 3 ) for 
a large a set of orbits selected from a uniform grid in initial condition space. Resonances 
appear on the frequency map as lines, either densely filled lines in the case of stable 
resonances, or gaps in the case of unstable resonances; the frequency map is effectively a 
representation of the Arnold web (Laskar 1993). Resonances are discussed in more detail 
in §3. 

2.3 Iterative Approaches 

Iterative approaches to torus construction consist of finding successively better approxi- 
mations to the map 9 — > x given some initial guess x(#); canonical perturbation theory is 
a special case, and in fact iterative schemes often reduce to perturbative methods in ap- 
propriate limits. Iterative algorithms were first developed in the context of semi-classical 




Figure 2. Construction of a 2 DOF, box-orbit torus in a Stdckel potential using the 
NAFF trajectory- following algorithm, (a) The orbit and its actions, computed using equa- 
tion (4) with k max terms. Dashed lines show the exact Jj. (b) The map y{0i,9 2 ); dashed 
contours correspond to negative values ofy. A(k max ) is the RMS error in the reconstructed 
map, calculated using an equation similar to (5). 



quantization for computing energy levels of bound molecular systems, and they are still 
best suited to assigning energies to actions, H(J). Most of the other quantities of interest 
to galactic dynamicists - e.g. the fundamental frequencies ooi - are not easily recovered 
using these algorithms. Iterative schemes also tend to be numerically unstable unless the 
initial guess is close to the true solution. On the other hand, iterative algorithms can 
be more efficient than trajectory- following algorithms for orbits that are near (but not 
exactly on) resonances. 

Ratcliff, Chang & Schwarzschild (1984) pioneered iterative schemes in galactic dynam- 
ics. They noted that the equations of motion of a 2 DOF regular orbit, 



can be written in the form 



d d V <9$ 



^W^^WJ x a*' 



If one specifies u\ and UJ2 and treats d^/dx and d$/dy as functions of the 8{, equations 
(0) can be viewed as nonlinear differential equations for x(6i, 82) and y(9i, 82). Ratcliff et 
al. expressed the coordinates as Fourier series in the angle variables, 



x(e) = ^X n e m -' 



Substituting into (|7|) gives 



n • uj) 2 X n e in - e = V$ (9) 



where the right hand side is again understood to be a function of the angles. Ratcliff et al. 
truncated the Fourier series after a finite number of terms and required equations to 
be satisfied on a grid of points around the torus. They then solved for the X n by iterating 
from an initial guess. Convergence was found to be possible if the initial guess was close 
to the exact solution. A similar algorithm was developed for recovering tori in the case 
that the actions, rather than the frequencies, are specified a priori. Guerra & Ratcliff 
(1990) applied these algorithms to motion in the plane of rotation of a nonaxisymmetric 
potential. 

Another iterative approach to torus construction was developed by Chapman, Garrett 
& Miller (1976) in the context of semiclassical quantum theory. One begins by dividing the 
Hamiltonian H into separable and non-separable parts H and Hi, then seeks a generating 
function S that maps the known tori of H into tori of H. For a generating function of 
the F 2 -type (Goldstein 1980), one has 

J(M') = §, ^> j/ ) = |f ( 10 ) 

where (J, 8) and (J', 8') are the action-angle variables of H and H respectively. The 
generator S is determined, for a specified J', by substituting the first of equations (|I0D 
into the Hamiltonian and requiring the result to be independent of 8. One then arrives 
at H(J'). Chapman et al. showed that a sufficiently general form for S is 

S(8,J')=8-J'- i J2S n (3')e me , (11) 

where the first term is the identity transformation, and they evaluated a number of it- 
erative schemes for finding the S n . One such scheme was found to recover the results 
of first-order perturbation theory after a single iteration. McGill & Binney (1990) ap- 
plied the Chapman et al. algorithm to 2 DOF motion in the axisymmetric logarithmic 
potential. 

The generating function approach is not naturally suited to deriving the other quan- 
tities of interest to galactic dynamicists. For instance, equation ( |10"|) gives 8' (8) as a 
derivative of S, but since S must be computed separately for every J' its derivative is 
likely to be ill-conditioned. Binney & Kumar (1993) and Kaasalainen & Binney (1994a) 
discussed two schemes for finding 8' (8); the first requires the solution of a formally infinite 



set of equations, while the latter requires multiple integrations of the equations of motion 
for each torus - effectively a trajectory-following scheme. 

Warnock (1991) presented a hybrid scheme in which the generating function S was 
derived by numerically integrating an orbit from appropriate initial conditions, transform- 
ing the coordinates to (J, 9) of Hq and interpolating J on a regular grid in 9. The values 
of the S n then follow from the first equation of ([H]) after a discrete Fourier transform. 
Kaasalainen & Binney (1994b) found that Warnock's scheme could be used to substan- 
tially refine the solutions found via their iterative algorithm. Another hybrid scheme was 
discussed by Reiman & Pomphrey (1991). 

Having computed the energy on a grid of J' values, one can interpolate to obtain the 
full Hamiltonian H(J'). If the system is not in fact completely integrable, this H may be 
rigorously interpreted as smooth approximation to the true H (Warnock & Ruth 1991, 
1992) and can be taken as the starting point for secular perturbation theory. Kaasalainen 
(1994) developed this idea and showed how to recover accurate surfaces of section in the 
neighborhood of low-order resonances in the planar logarithmic potential. 

Percival (1977) described a variational principle for constructing tori. His technique 
has apparently not yet been implemented in the context of galactic dynamics. 



2.4 Chaotic Motion 

Torus-construction machinery may be applied to orbits that are approximately, but not 
precisely, regular (Laskar 1993). The frequency spectrum of a weakly chaotic orbit will 
typically be close to that of a regular orbit, with most of the lines well approximated 
as linear combinations of three "fundamental frequencies" Ui. However these frequencies 
will change with time as the orbit migrates from one "torus" to another. The diffusion 
rate can be measured via quantities like 5uj = \ui — uj[\, the change in a "fundamental 
frequency" over two consecutive integration intervals. Papaphilippou & Laskar (1996, 
1998), Valluri & Merritt (1998) and Wachlin & Ferraz-Mello (1998) used this technique 
to study chaos and diffusion in triaxial galactic potentials. 

Measuring chaos via quantities like Su has a number of advantages over the traditional 
technique based on computation of the Liapunov exponents (Lichtenberg & Lieberman 
1992). 5u can be accurately determined after just a few tens of orbital periods, whereas 
determination of the Liapunov exponents may require much longer integrations. The 
Liapunov exponents measure only the rate of growth of infinitesimal perturbations around 
the trajectory, while 5u measures the finite "movement" of the trajectory in action-angle 
space, a more physically interesting measure of chaos. It is possible for orbits to be 
extremely unstable in the sense of having large Liapunov exponents, but to behave nearly 
regularly in the sense of having small 5uj\ an example is presented in §6. 



3 RESONANCES 



The character of a regular orbit depends critically on whether the fundamental frequencies 
Ui are independent, or whether they satisfy one or more nontrivial linear relations of the 



form 

JV 

5^771^ = (12) 

i=l 

with N the number of degrees of freedom and m ; integers, not all of which are zero. 
Generally there exists no relation like equation (|12|)> the frequencies are incommensurate, 
and the trajectory fills its invariant torus uniformly and densely in a time-averaged sense. 
When one or more resonance relations are satisfied, however, the trajectory is restricted 
to a phase-space region of lower dimensionality than N. 




Figure 3. Resonant tori, (a) A two-dimensional torus, shown here as a square with 
identified edges. The plotted trajectory satisfies a 2 : 1 resonance between the fundamental 
frequencies, Ui — 2u 2 = (e.g. a "banana"), (b) A three-dimensional torus, shown here as 
a cube with identified sides. The shaded region is covered densely by a resonant trajectory 
for which 2u)\ + ui 2 — 2uo^ = 0. This trajectory is not closed, but it is restricted by the 
resonance condition to a two-dimensional subset of the torus. The orbit in configuration 
space is thin. 

In the case of a two-dimensional regular orbit, the angle variables are 

0i = wit + 010, e 2 = u 2 t + e 20 , (13) 

which define the surface of a torus (Fig. 1). Because of the quasi-periodicity of the orbit, 
its torus can be mapped onto a square in the (0i, 02)-plane, with each side ranging from 
to 27r (Figure 3a); the top and bottom of the square are identified with each other, as are 
the left and right sides. In the general case, the frequencies u\ and u 2 are incommensurate 
and the trajectory densely covers the entire (0i, 02)-plane after an infinite time. However 
if the ratio uj\/uj 2 = \m 2 jm\\ is a rational number, i.e. if m\ and m 2 are integers, the 
orbit closes on itself after \m 2 \ revolutions in 0i and \m\\ revolutions in 6 2 and fills only 
a one-dimensional subset of its torus (e.g. Arnold 1963, p. 164). Its dimensionality in 
configuration space is also one - the orbit is closed. Such an orbit has a single fundamental 
frequency ojq = Wi/m 2 = io 2 jm\ = 2n/T, with T the orbital period; after an elapsed time 
T, the trajectory returns to its starting point in phase space. Examples of resonant orbits 
in two-dimensional galactic potentials are the "boxlets" (Miralda-Escude & Schwarzschild 
1989). 



In the case of a three-dimensional regular orbit, the angle variables are 



X = (Jit + 0io, 02 = UJ 2 t + 6 20 , 03 = U 3 t + 3o . (14) 

The orbit may now be mapped into a cube whose axes are identified with the 0i (Figure 
3b). If the Ui are incommensurate, this cube will be densely filled after a long time. 
However if a single condition of the form 

irtiuj 1 + m 2 uj2 + m 3 uj 3 = (15) 

is satisfied with integer rrii, the motion is restricted for all time to a two-dimensional 
subset of its torus . Such an orbit is not closed; instead, as suggested by Figure 3b, it is 
thin, confined to a sheet or membrane in configuration space, which it fills densely after 
infinite time. 



Just as in the two-dimensional case, the condition ( jT5|) may be used to reduce the 
number of independent frequencies by one. Defining the two "base" frequencies 0Uq , Uq 

as 

u<p = u) 3 /m 1 , uf> = LU 2 /mi, (16) 



we may write 

(1) (2) 

(2) 

u 3 = miu^\ (17) 
Since the motion is quasi-periodic, i.e. 

x(t) = 2J x fc ex P i {h^i + m k u} 2 + n k u 3 ) t, (18) 

k 

with (4, rrik, n^.) integers, it will remain quasi-periodic when expressed in terms of the two 
base frequencies: 

x(t) = ^X fc expi (-l k m 3 + n k m l )uj^ ) + {-l k m 2 + m fc mi) uj^ ] t 

k 

= ^ X fc exp i (ik'uj^ + rrik'ujQ^ t, 
k 

= ^Xfeexp^/^W+m^ 2 )), 
k 

Ik = -hm 3 + n k mi, m k ' = -l k m 2 + m k mi, 
0W = 4% 0^=J*h. (19) 

A Fourier transform of the motion will therefore consist of a set of spikes whose locations 
can be expressed as linear combinations of just two frequencies. Equation (|T9| ) is a para- 
metric expression for the Cartesian coordinates in terms of the angles on the 2-torus, i.e. 
it is a reconstruction of the (reduced) torus. A number of examples of resonant box orbits 
reconstructed in this way are illustrated in Figures 4 and 5. 



Certain special orbits may satisfy two independent resonance relations simultaneously. 
In this case we can write: 

m-iUi + m 2 uj 2 + "^3^3 = 0, 
niui + n 2 uo 2 + n 3 uo 3 = 0, 

and each frequency u>i may be expressed as a rational fraction of any other: 

loi _ m 2 n 3 - m 3 n 2 _ h _ m 3 n x - m x n 3 _ k 

uo 3 m\n 2 - m 2 ni l 3 u 3 mxn 2 -m 2 nx l 3 

with (h,l 2 ,l 3 ) integers. The motion is therefore periodic with a single base frequency 
u —0Ji/li = uj 2 jl 2 = uj 3 /1 3 and the trajectory is closed - the orbit is a three-dimensional, 
closed curve. In a system with N degrees of freedom, N — 1 such conditions are required 
for closure; only in the 2DOF case does a single resonance condition imply closure. 

Following Poincare (1892), it has commonly been assumed that closed orbits are the 
fundamental "building blocks" of phase space. However in three-dimensional potentials, 
one expects thin orbits to be more common than closed ones, in the sense that orbits 
satisfying one resonance condition are more likely than orbits satisfying two. Hence one 
expects that most regular orbits will be associated with families whose parent is a thin 
orbit. Numerical integrations of orbits in realistic nonaxisymmetric potentials suggest 
that this is in fact the case: the majority of regular orbits have most of their "power" 
in frequencies that lie close to linear combinations of two fundamental frequencies (thin 
orbit) rather than one frequency (closed orbit) (Merritt & Valluri 1999; Figure 7). 



(20) 



4 ORBITAL STRUCTURE OF TRIAXIAL POTEN- 
TIALS WITH CENTRAL SINGULARITIES 

Non-integrability is likely to be a generic feature of galactic potentials, for two reasons. 
First, galaxies are often observed to be non-axisymmetric, either due to the presence 
of imbedded sub-systems like bars, or because the stellar distribution is globally triax- 
ial. Observational evidence for global triaxiality in elliptical galaxies is not particularly 
strong; few ellipticals exhibit significant minor-axis rotation (Franx, Illingworth & de 
Zeeuw 1991), and detailed modelling of a handful of nearby ellipticals suggests that their 
kinematics can often be very well reproduced by assuming axisymmetry (e.g. van der 
Marel et al. 1998). However, at least some elliptical galaxies and bulges exhibit clear 
kinematical signatures of non-axisymmetry (s.g. Schechter & Gunn 1979; Franx, Illing- 
worth & Heckman 1989), and the observed distribution of Hubble types is likewise incon- 
sistent with the assumption that all ellipticals are precisely axisymmetric (Tremblay & 
Merritt 1995, 1996; Ryden 1996). Mergers between disk galaxies also produce generically 
triaxial systems (Barnes 1996), and departures from axisymmetry (possibly transient) 
are widely argued to be necessary for the rapid growth of nuclear black holes during the 
quasar epoch (Shlosman, Begelman & Frank 1990), for the fueling of starburst galaxies 
(Sanders & Mirabel 1996), and for the large radio luminosities of some ellipticals (Bicknell 
et al. 1997). These arguments suggest that most elliptical galaxies or bulges may have 
been triaxial at an earlier epoch, and perhaps that triaxiality is a recurrent phenomenon 
induced by mergers or other interactions. 



The second feature of galactic potentials that is conducive to non-integrability is the 
apparently universal presence at the centers of stellar spheroids of high stellar densities 
and supermassive black holes. Low-luminosity ellipticals and bulges have stellar lumi- 
nosity profiles that diverge as unbroken power laws at small radii, p ~ r -7 , with 7^2. 
Brighter galaxies also exhibit power laws in the space density of stars, but with shallower 
slopes, 7 < 1; seen in projection, these weaker cusps appear as cores (Kormendy 1985). 
The gravitational force in an r~ 2 density cusp diverges as r _1 , not steep enough to produce 
large-angle deflections in the motion of stars that pass near the center. However galaxies 
also contain supermassive black holes, with masses that correlate astonishingly well with 
the velocity dispersion of the stars (Ferrarese & Merritt 2000); the ratio of black hole 
mass to spheroid mass is ~ 0.0015 with small scatter (Merritt & Ferrarese 2001). The 
combination of non-axisymmetry in the potential with a steep central force gradient is 
conducive to non-integrability and chaos, since many orbits in non-axisymmetric poten- 
tials pass near the center where they undergo strong gravitational deflections (Gerhard & 
Binney 1985). 

In a triaxial potential containing a central point mass, the phase space divides naturally 
into three regions depending on energy, i.e. on distance from the center (Figure 6). In 
the innermost region, where the enclosed mass in stars is less than the mass of the black 
hole, the potential is dominated by the central singularity and the motion is essentially 
regular. The gravitaional force from the stars acts as a small perturbation causing the 
nearly-Keplerian orbits around the black hole to slowly precess. The two major orbit 
families in this region are the tube orbits, high angular momentum orbits that avoid the 
center; and the pyramid orbits, Keplerian ellipses that precess in two orthogonal planes 
parallel to the short axis of the figure (Sridhar & Touma 1999; Sambhus & Sridhar 2000; 
Poon & Merritt 2001). Pyramid orbits are similar to the classical box orbits of integrable 
triaxial potentials except that their elongation is counter to that of the triaxial figure, 
making them less useful for self-consistently reconstructing a galaxy's shape. 

At intermediate radii, the black hole acts as a scattering center rendering almost all 
of the center-filling or box orbits stochastic. (Tube orbits persist at these and higher 
energies and remain mostly regular.) This "zone of chaos" extends from a few times r g , 
the radius where the black hole dominates the gravitational force, out to a radius where 
the enclosed stellar mass is roughly 10 2 times the mass of the black hole. The transition 
to chaos at r > r g is very rapid and occurs at lower energies in more elongated potentials 
(Poon & Merritt 2001). 

If the black hole mass exceeds ~ 10~ 2 times the mass of the stellar spheroid, as it may 
do in a few galaxies (Merritt & Ferrarese 2001), the chaotic zone will include essentially 
the entire potential outside of ~ r g . However if M, ps 10~ 3 M gal , as in the majority of 
galaxies, there exists a third, outermost region where the phase space is a complex mixture 
of chaotic and regular trajectories, including resonant box orbits like those in Figures 4 and 
5 that remain stable by avoiding the center (Carpintero & Aguilar 1998; Papaphillipou & 
Laskar 1998; Valluri & Merritt 1998; Wachlin & Ferraz-Mello 1998). Figure 7 illustrates 
the complexity of box-orbit phase space at large energies in two triaxial potentials: one 
with a weak density cusp and the other with a central point mass. 

Non-integrable potentials often exhibit a transition to global stochasticity as the mag- 
nitude of some perturbation parameter is increased. The results summarized above sug- 
gest that there are two such perturbation parameters associated with motion in triaxial 



galaxies containing central black holes. In a triaxial galaxy with a given M., the motion 
of center-filling orbits undergoes a sudden transition to stochasticity as the energy is in- 
creased; the critical value is the energy at which the gravitational force from the stars is of 
order the force from the black hole. If one imagines increasing M. in an otherwise fixed, 
triaxial potential, the zone of chaos that extends outward from this radius will eventually 
encompass the entire potential; this occurs when the second "perturbation parameter," 
M./Mg a i, exceeds ~ 1CT 2 . Thus at intermediate radii, in the "zone of chaos," and per- 
haps throughout an elliptical galaxy containing a central black hole, triaxiality should be 
difficult to maintain. 



5 MIXING AND COLLISIONLESS RELAXATION 



Stochastic motion introduces a new time scale into galactic dynamics, the mixing time. 
Mixing is the process by which a non-uniform distribution of particles in phase space 
relaxes to a uniform distribution, at least in a coarse-grained sense. A weak sort of mixing, 
phase mixing, occurs even in integrable potentials, as particles on adjacent tori gradually 
move apart (Lynden-Bell 1967; Figure 8a). Phase mixing is responsible for the fact that 
the coarse-grained phase space density in relaxed integrable systems is nearly constant 
around tori. A stronger sort of mixing takes place in chaotic systems. Chaotic motion is 
essentially random in the sense that the likelihood of finding a particle anywhere in the 
stochastic region tends toward a constant value after a sufficiently long time. An initially 
compact group of stars should therefore spread out until it covers the accessible phase 
space region uniformly in a coarse-grained sense (Kandrup & Mahon 1994; Figure 8b). 
This "chaotic mixing" is irreversible in the sense that an infinitely fine tuning of velocities 
would be required in order to undo its effects. It also occurs on a characteristic time 
scale, the Liapunov time associated with exponential divergence of nearby trajectories. 
Phase mixing, by contrast, has no associated time scale; its rate depends on the range of 
frequencies associated with orbits in the region of interest, and this rate tends to zero in 
the case of a set of trajectories drawn from a single invariant torus - a set of points on 
the torus translates, unchanged, around the torus. 

Figure 9 shows examples of chaotic mixing in a triaxial potential with a central point 
mass. Ensembles of orbits were started at rest on an equipotential surface and integrated 
in tandem for several crossing times. The central point had a mass M. = 0.03 in units of 
the galaxy mass. The first ensemble (a) was begun on an equipotential surface enclosing 
a mass ~ 3 times that of the central point; for ensembles (b) and (c) these ratios were 
~ 7 and ~ 17 respectively - all within the "zone of chaos." Mixing occurs very rapidly 
in these ensembles. At the lowest energy (ensemble a), the linear extent of the points in 
configuration space roughly doubles every crossing time until T ss 4, when the volume 
defined by the equipotential surface appears to be nearly filled. At the highest energy (en- 
semble c), mixing is slower but substantial changes still take place in a few crossing times. 
The final distribution of points at this energy still shows some structure, reminiscent of a 
box orbit. 

The irreversibility of mixing flows like the ones illustrated in Figure 9 implies a re- 
duction in the effective number of orbits: all the stochastic trajectories at a given energy 
are gradually replaced by a single invariant ensemble, whose shape is typically not well 



matched to that of the galaxy (Merritt & Fridman 1996). If time scales for chaotic mix- 
ing are comparable to galaxy lifetimes, this reduction might be expected to encourage a 
galaxy to evolve away from a triaxial shape toward a more axisymmetric one, in which 
most of the orbits are tubes that avoid the destabilizing center. Such evolution has in fact 
been observed in iV-body simulations of the response of a triaxial galaxy to the growth 
of a central black hole. Merritt & Quinlan (1998) found that a triaxial galaxy evolves to 
axisymmetry in little more than the local crossing time at each radius when the black 
hole mass exceeds ~ 2.5% of the total galaxy mass. This is about an order of magnitude 
larger than the typical black hole mass ratio in real galaxies (Merritt & Ferrarese 2001), 
but Merritt & Quinlan observed more gradual evolution even when the mass ratio was 10 
times smaller, at a rate that would imply substantial shape changes over a galaxy lifetime. 
These simulations suggest an explanation for the generally low level of triaxiality observed 
in real galaxies (Bak & Statler 2000). 



6 CHAOS IN COLLISIONAL SYSTEMS 

The discussion presented so far has assumed that galaxy potentials are smooth, or "col- 
lisionless." In fact, the gravitational force on a star in a galaxy can be broken up into 
two components: a rapidly varying component that arises from the discrete distribu- 
tion of stars, and a smoothly varying component that arises from the large-scale matter 
distribution. The effects of the discrete component of the force relative to the smooth 
component are usually assumed to scale as ~ In N/N, the ratio of dynamical to two-body 
relaxation times. For galaxies, which have N ~ 10 11 , collisional effects should therefore 
be unimportant. 

If this were the case, it should be possible to show that the iV-body trajectories go 
over, in the limit of large N, to the orbits in the corresponding smoothed-out potential - 
i.e., that the equations of motion of the iV-body problem tend to the characteristics of the 
collisionless Boltzmann equation as N — > oo. However this has never been demonstrated, 
and in fact there is an important sense in which the equations of motion in an iV-body 
system do not tend toward the trajectories of the corresponding smooth potential in the 
limit of large N. 

This surprising statement is justified in Figure 10, which shows the results of test- 
particle integrations in a potential consisting of iV fixed point masses distributed randomly 
and uniformly within a triaxial ellipsoid. The mass of each of the N points is m = 1/N, 
so that the total mass and mean density of the ellipsoid remain constant as N is varied. 
In the limit N — > oo, one might expect the equations of motion to approach those of 
a 3d harmonic oscillator, since the potential of a uniform ellipsoid is quadratic in the 
coordinates, $ = $o^Aii (Chandrasekhar 1969). However the upper left-hand panel 
shows that the Liapunov exponents a of orbits in the iV-body potential do not tend to 
zero with increasing N. Instead, the instability time scale appears to r each a roughly 
constant value (expressed as a fraction of the crossing time T cr ) for N > 10 3 . Furthermore 
the instability time scale is very short, a fraction of the crossing time! 

The generic instability of the iV-body problem was first noted by Miller (1964), who 
calculated the time evolution of the separation between two iV-body systems with slightly 



different initial conditions. He defined this separation as 

A(f) = [£ (x 2 - Xl ) 2 + ^ (v 2 - Vl ) 2 ] V2 

with xi and x 2 the N configuration-space coordinates in iV-body systems 1 and 2 and 
vi and v 2 the velocities; the summations extend over all the particles. Miller found, for 
4 < iV < 32, that A grew roughly exponentially with a characteristic time scale that was 
a fraction of the crossing time, as in the fixed iV-body problem of Figure 10. 

What are the physical implications of this generic instability? Several suggestions have 
been made. Gurzadyan & Saviddy (1986), who first investigated the large- N dependence 
of the instability using an idealized model, suggested that the exponential divergence 
implies chaotic mixing on a similar time scale, and hence that stellar systems should relax 
much more rapidly than implied by the standard Chandrasekhar formula. Heggie (1991) 
disagreed, but suggested that the use of smooth potentials for approximating galaxies 
would need to be abandoned, at least for studies of orbital instability. Kandrup (1998) 
suggested that - while individual orbits may always be exponentially unstable - ensembles 
of iV-body systems might behave, on average, as if the potential were smooth. 

Figure 10 suggests an even stronger way in which the motion goes over to that of 
the collisionless problem as N — > oo. The open circles in the upper left-hand panel of 
that figure show a second measure of the orbital evolution: the rms variation, over 20 
orbital periods, of the action J x for each ensemble of orbits. Contrary to the behavior 
of the Liapunov exponents, the average changes in the actions tend uniformly to zero 
as N is increased - in other words, the orbits approach more and more closely, in their 
macroscopic behavior, to that of integrable orbits even though they remain locally unstable 
(as measured by the Liapunov exponents) to a degree that is nearly independent of N. 
Plots of the trajectories of some typical orbits (lower left panel) confirm this interpretation. 
These results suggest the way in which trajectories in the iV-body problem tend toward 
those in the corresponding smooth potential: as iV is increased, orbits are confined more 
and more strongly to narrow regions of phase space around the invariant tori of the smooth 
potential. 

It is remarkable that orbits can be extremely unstable locally, as measured by their 
Liapunov exponents, and yet behave macroscopically in a way that is essentially identical 
to that of regular orbits. Apparently, the exponential growth of perturbations must 
saturate at some finite amplitude, and this saturation amplitude must be a decreasing 
function of N. The lower right-hand panel of Figure 10 verifies this conjecture for a 
few pairs of orbits with nearly identical initial conditions. The early divergence takes 
place at a rate that is independent of N, but for large N, the separation saturates at a 
value that is much smaller than the size of the system. These pairs of orbits act as if 
they are confined to the same, restricted region of phase space; saturation occurs when the 
separation between them is of order the width of this region. The fact that the exponential 
divergence saturates sooner for larger N suggests that the width of the confining regions 
decreases with increasing N. 

These results suggest that collisional relaxation in stellar systems is intimately con- 
nected with the evolution of orbits under conditions of weak chaos, i.e., with Arnold 
diffusion. This connection would be a fruitful topic for future study. 

Some of the work presented here was first published in collaboration with M. Valluri. 
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Figure 4. Surfaces filled by a set of thin, or resonant, box orbits in a non-integrable 
triaxial potential (Merritt & Valluri 1999), as seen from vantage points on each of the 
three principal axes. The cross sections of these orbits are shown in Figure 5. 
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Figure 5. Intersections with the principal planes of the thin box orbits shown in Figure 
4- Because the orbits are thin, their intersections with any plane define a curve or set of 
curves. The center of the potential is indicated by a cross. 




Figure 6. Three zones in the phase space of triaxial potentials (see text). 
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Figure 7. Non-integrability in triaxial potentials (Merritt & Valluri 1999). The mass 
model in (a) has a weak (7 = 0.5 ^ density cusp and no black hole; in (b) the black hole 
contains 0.3% of the total mass. Each panel shows one octant of an equipotential surface, 
lying close to the half-mass radius of the model; the z (short) axis is vertical and the x 
(long) axis is to the left. The grey scale measures the degree of stochasticity of orbits 
started with zero velocity on the equipotential surface. Stable resonance zones - the white 
bands in (a) and (b) - are labelled by their defining integers (m 1 ,m 2 ,m 3 ). Panels (c) and 
(d) show the pericenter distance A of a set of 10 3 orbits with starting points along the 
heavy solid lines in (a) and (b). Panels (e) and (f) plot a measure of the chaos for these 
orbits; 5u/u is the fractional change in the frequency of the strongest line in the orbit's 
frequency spectrum. 




Figure 8. (a) Phase mixing vs. (b) chaotic mixing. 



Figure 9. Mixing in a triaxial potential with a central point containing 3% of the total 
mass (Valluri & Merritt 2000). Time is in units of the local crossing time. Ensembles of 
10 4 phase points were distributed initially (T = 0) in patches on a equipotential surface 
with zero velocity. 
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Figure 10. Evolution of orbits in a potential consisting of N fixed point masses, 
m = 1/N, distributed randomly and uniformly in an ellipsoidal volume (Valluri & Merritt 
2000) (see text). 



